### Daniel J. Hopkins
### Replication code:
### "The Activation of Prejudice and Presidential Voting:
### Panel Evidence from the 2016 U.S. Election"
### Conditionally accepted, Political Behavior
### July 31, 2019

#### load libraries ----
library(texreg)
library(xtable)
library(MNP)
library(lme4)
library(pBrackets)

### useful functions
is.zero <- function(x){
	1*(x==0)
}
is.one <- function(x){
	1*(x==1)
}
is.two <- function(x){
	1*(x==2)
}


#### load original data ----
#### note: the original data is not provided
#### but the code below is provided to allow 
#### users to understand variable construction
# load("2016wave/merged-2014-2016-02012018.Rdata")
# 
# ##### recode variables ----
# dta.all3$WBPREJUDICE45 <- apply(cbind(dta.all3$WBPREJUDICE4,dta.all3$WBPREJUDICE5),1,mean)/100
# 
# dta.all3$WBPREJUDICE4R <- dta.all3$WBPREJUDICE4/100
# dta.all3$WBPREJUDICE5R <- dta.all3$WBPREJUDICE5/100
# 
# dta.all3$WBPREJUDICE6R <- dta.all3$WBPREJUDICE6/100
# dta.all3$WBPREJUDICE7R <- dta.all3$WBPREJUDICE7/100
# dta.all3$WBPREJUDICE10R <- dta.all3$WBPREJUDICE10/100
# dta.all3$WBPREJUDICE11R <- dta.all3$WBPREJUDICE11/100
# dta.all3$WBPREJUDICE12R <- dta.all3$WBPREJUDICE12/100
# 
# dta.all3$WHPREJUDICE6R <- dta.all3$WHPREJUDICE6/100
# dta.all3$WHPREJUDICE7R <- dta.all3$WHPREJUDICE7/100
# dta.all3$WHPREJUDICE10R <- dta.all3$WHPREJUDICE10/100
# dta.all3$WHPREJUDICE11R <- dta.all3$WHPREJUDICE11/100
# 
# dta.all3$WBPREJUDICE100 <- dta.all3$WBPREJUDICE10
# dta.all3$WBPREJUDICE100[dta.all3$STTRSBLK100==1 & ! dta.all3$STTRSBLK100 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE100[dta.all3$STHRDBLK100==1 & ! dta.all3$STHRDBLK100 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE100[dta.all3$STTRSWHT100==1 & ! dta.all3$STTRSWHT100 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE100[dta.all3$STHRDWHT100==1 & ! dta.all3$STHRDWHT100 %in% c(NA)] <- NA
# 
# dta.all3$WBPREJUDICE110 <- dta.all3$WBPREJUDICE11
# dta.all3$WBPREJUDICE110[dta.all3$STTRSBLK110==1 & ! dta.all3$STTRSBLK110 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE110[dta.all3$STHRDBLK110==1 & ! dta.all3$STHRDBLK110 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE110[dta.all3$STTRSWHT110==1 & ! dta.all3$STTRSWHT110 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE110[dta.all3$STHRDWHT110==1 & ! dta.all3$STHRDWHT110 %in% c(NA)] <- NA
# 
# dta.all3$WBPREJUDICE120 <- dta.all3$WBPREJUDICE12
# dta.all3$WBPREJUDICE120[dta.all3$STTRSBLK120==1 & ! dta.all3$STTRSBLK120 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE120[dta.all3$STHRDBLK120==1 & ! dta.all3$STHRDBLK120 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE120[dta.all3$STTRSWHT120==1 & ! dta.all3$STTRSWHT120 %in% c(NA)] <- NA
# dta.all3$WBPREJUDICE120[dta.all3$STHRDWHT120==1 & ! dta.all3$STHRDWHT120 %in% c(NA)] <- NA
# 
# dta.all3$WHPREJUDICE100 <- dta.all3$WHPREJUDICE10
# dta.all3$WHPREJUDICE100[dta.all3$STTRSHIS100==1 & ! dta.all3$STTRSHIS100 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE100[dta.all3$STHRDHIS100==1 & ! dta.all3$STHRDHIS100 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE100[dta.all3$STTRSWHT100==1 & ! dta.all3$STTRSWHT100 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE100[dta.all3$STHRDWHT100==1 & ! dta.all3$STHRDWHT100 %in% c(NA)] <- NA
# 
# dta.all3$WHPREJUDICE110 <- dta.all3$WHPREJUDICE11
# dta.all3$WHPREJUDICE110[dta.all3$STTRSHIS110==1 & ! dta.all3$STTRSHIS110 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE110[dta.all3$STHRDHIS110==1 & ! dta.all3$STHRDHIS110 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE110[dta.all3$STTRSWHT110==1 & ! dta.all3$STTRSWHT110 %in% c(NA)] <- NA
# dta.all3$WHPREJUDICE110[dta.all3$STHRDWHT110==1 & ! dta.all3$STHRDWHT110 %in% c(NA)] <- NA
# 
# dta.all <- dta.all3
# 
# dta.all$STRBLACK6 <- apply(cbind(dta.all$STHRDBLK6,dta.all$STINTBLK6,dta.all$STTRSBLK6),1,mean)/100
# dta.all$STRBLACK7 <- apply(cbind(dta.all$STHRDBLK7,dta.all$STINTBLK7,dta.all$STTRSBLK7),1,mean)/100
# 
# dta.all$STRBLACK67 <- apply(cbind(dta.all$STRBLACK6,dta.all$STRBLACK7),1,mean)
# 
# dta.all$WBPREJUDICE67 <- apply(cbind(dta.all$WBPREJUDICE6,dta.all$WBPREJUDICE7),1,mean)/100
# dta.all$WHPREJUDICE67 <- apply(cbind(dta.all$WHPREJUDICE6,dta.all$WHPREJUDICE7),1,mean)/100
# 
# ##
# dta.all$TRUMPVC10F <- as.factor(dta.all$TRUMPVC10B)
# dta.all$TRUMPVC11F <- as.factor(dta.all$TRUMPVC11B)
# dta.all$TRUMPVC12F <- as.factor(dta.all$TRUMPVC12B)
# 
# dta.all$ROMNEYVO6F <- as.factor(dta.all$ROMNEYVO6B)
# dta.all$ROMNEYVO7F <- as.factor(dta.all$ROMNEYVO7B)
# 
# dta.all$MCCAINVO4A1 <- 1*(dta.all$MCCAINVO4A==1)
# dta.all$MCCAINVO4A9 <- 1*(dta.all$MCCAINVO4A==9)
# 
# dta.all$ROMNEYVO6A1 <- 1*(dta.all$ROMNEYVO6A==1)
# dta.all$ROMNEYVO6A9 <- 1*(dta.all$ROMNEYVO6A==9)
# 
# dta.all$WEAKDEM0 <- 1*(dta.all$PID0==2)
# dta.all$INDDEM0 <- 1*(dta.all$PID0==3)
# dta.all$INDEP0 <- 1*(dta.all$PID0==4)
# dta.all$INDGOP0 <- 1*(dta.all$PID0==5)
# dta.all$WEAKGOP0 <- 1*(dta.all$PID0==6)
# dta.all$STRONGGOP0 <- 1*(dta.all$PID0==7)
# 
# dta.all$WEAKDEM4 <- 1*(dta.all$PID4==2)
# dta.all$INDDEM4 <- 1*(dta.all$PID4==3)
# dta.all$INDEP4 <- 1*(dta.all$PID4==4)
# dta.all$INDGOP4 <- 1*(dta.all$PID4==5)
# dta.all$WEAKGOP4 <- 1*(dta.all$PID4==6)
# dta.all$STRONGGOP4 <- 1*(dta.all$PID4==7)
# 
# dta.all$WEAKDEM6 <- 1*(dta.all$PID6==2)
# dta.all$INDDEM6 <- 1*(dta.all$PID6==3)
# dta.all$INDEP6 <- 1*(dta.all$PID6==4)
# dta.all$INDGOP6 <- 1*(dta.all$PID6==5)
# dta.all$WEAKGOP6 <- 1*(dta.all$PID6==6)
# dta.all$STRONGGOP6 <- 1*(dta.all$PID6==7)
# 
# dta.all$WEAKDEM10 <- 1*(dta.all$PID10==2)
# dta.all$INDDEM10 <- 1*(dta.all$PID10==3)
# dta.all$INDEP10 <- 1*(dta.all$PID10==4)
# dta.all$INDGOP10 <- 1*(dta.all$PID10==5)
# dta.all$WEAKGOP10 <- 1*(dta.all$PID10==6)
# dta.all$STRONGGOP10 <- 1*(dta.all$PID10==7)
# 
# dta.all$BUSHVK0B <- NA
# dta.all$BUSHVK0B[dta.all$pppa0242_1=="George W. Bush (Republican)"] <- 1
# dta.all$BUSHVK0B[dta.all$pppa0242_1=="John Kerry (Democrat)"] <- -1
# dta.all$BUSHVK0B[dta.all$pppa0242_1 %in% c("Ralph Nader (Green)","Another candidate, please specify:")] <- 0
# dta.all$BUSHVK0B[dta.all$pppa0241_1 %in% c("No")] <- 0
# 
# dta.all$BUSHVK0A1 <- 1*(dta.all$BUSHVK0B==1)
# dta.all$BUSHVK0A9 <- 1*(dta.all$BUSHVK0B==0)
# 
# dta.all$GOPVDEM06B <- NA
# dta.all$GOPVDEM06B[dta.all$pppa0302_1=="The Republican candidate"] <- 1
# dta.all$GOPVDEM06B[dta.all$pppa0302_1=="The Democratic candidate"] <- -1
# dta.all$GOPVDEM06B[dta.all$pppa0302_1 %in% c("A candidate from another party")] <- 0
# dta.all$GOPVDEM06B[dta.all$pppa0301_1 %in% c("No")] <- 0
# 
# dta.all$MCCAINVO3B <- NA
# dta.all$MCCAINVO3B[dta.all$MCCAINVO3A==0] <- -1
# dta.all$MCCAINVO3B[dta.all$MCCAINVO3A==9] <- 0
# dta.all$MCCAINVO3B[dta.all$MCCAINVO3A==1] <- 1
# 
# ### SUMMARIZE DATA
# 
# dta.all$HASHS <- 1*(dta.all$EDYEARS6 >= 12)
# dta.all$HASBA <- 1*(dta.all$EDYEARS6 >= 16)
# dta.all$OVER65 <- 1*(dta.all$AGE6 >= 65)
# 
# dta.all$WBESTEEM67 <- dta.all$WBPREJUDICE67
# dta.all$WBESTEEM67[dta.all$WBPREJUDICE67 > 0] <- 0
# 
# dta.all$WBPREJUDICEONLY67 <- dta.all$WBPREJUDICE67
# dta.all$WBPREJUDICEONLY67[dta.all$WBPREJUDICE67 < 0] <- 0
# 
# ##### issue attitudes for factor analysis
# 
# dta.all$BUILDWALL4R <- dta.all$BUILDWALL1R <- NA
# dta.all$BUILDWALL4R[dta.all$BUILDWALL4==1] <- 4
# dta.all$BUILDWALL4R[dta.all$BUILDWALL4==2] <- 3
# dta.all$BUILDWALL4R[dta.all$BUILDWALL4==3] <- 2
# dta.all$BUILDWALL4R[dta.all$BUILDWALL4==4] <- 1
# 
# dta.all$BUILDWALL1R[dta.all$BUILDWALL1==1] <- 4
# dta.all$BUILDWALL1R[dta.all$BUILDWALL1==2] <- 3
# dta.all$BUILDWALL1R[dta.all$BUILDWALL1==3] <- 2
# dta.all$BUILDWALL1R[dta.all$BUILDWALL1==4] <- 1
# 
# dta.all$IMMIGIDX <- (dta.all$PATHWAY6-mean(dta.all$PATHWAY6,na.rm=T)/sd(dta.all$PATHWAY6,na.rm=T))+(dta.all$PATHWAY7-mean(dta.all$PATHWAY7,na.rm=T)/sd(dta.all$PATHWAY7,na.rm=T))+(dta.all$BUILDWALL1-mean(dta.all$BUILDWALL1R,na.rm=T)/sd(dta.all$BUILDWALL1R,na.rm=T))+(dta.all$BUILDWALL4R-mean(dta.all$BUILDWALL4R,na.rm=T)/sd(dta.all$BUILDWALL4R,na.rm=T))
# 
# dta.all$NAFTA7 <- NA
# dta.all$NAFTA7[dta.all$POS7_7=="Strongly favor"] <- 4
# dta.all$NAFTA7[dta.all$POS7_7=="Somewhat favor"] <- 3
# dta.all$NAFTA7[dta.all$POS7_7=="Somewhat oppose"] <- 2
# dta.all$NAFTA7[dta.all$POS7_7=="Strongly oppose"] <- 1
# 
# dta.all$NAFTA6 <- NA
# dta.all$NAFTA6[dta.all$POS7_6=="Strongly favor"] <- 4
# dta.all$NAFTA6[dta.all$POS7_6=="Somewhat favor"] <- 3
# dta.all$NAFTA6[dta.all$POS7_6=="Somewhat oppose"] <- 2
# dta.all$NAFTA6[dta.all$POS7_6=="Strongly oppose"] <- 1
# 
# dta.all$NAFTA4R <- NA
# dta.all$NAFTA4R[dta.all$NAFTA4==1] <- 4
# dta.all$NAFTA4R[dta.all$NAFTA4==2] <- 3
# dta.all$NAFTA4R[dta.all$NAFTA4==3] <- 2
# dta.all$NAFTA4R[dta.all$NAFTA4==4] <- 1
# 
# dta.all$GAYMAR7 <- NA
# dta.all$GAYMAR7[dta.all$POS1_8=="I support full marriage rights for gay and lesbian couples."] <- 3
# dta.all$GAYMAR7[dta.all$POS1_8=="I support civil unions or domestic partnerships, but not gay marriage."] <- 2
# dta.all$GAYMAR7[dta.all$POS1_8=="I do not support any form of legal recognition of the relationships of gay and lesbian couples."] <- 1
# 
# dta.all$PROCHOICE4 <- NA
# dta.all$PROCHOICE4[dta.all$POS4_4=="Abortion should be available to anyone who wants it."] <- 4
# dta.all$PROCHOICE4[dta.all$POS4_4=="Abortion should be available, but with stricter limits than it is now."] <- 3
# dta.all$PROCHOICE4[dta.all$POS4_4=="Abortion should not be permitted except in cases of rape, incest, or when the life of the woman is at risk."] <- 2
# dta.all$PROCHOICE4[dta.all$POS4_4=="Abortion should not be permitted under any circumstances."] <- 1
# 
# dta.all$PROCHOICE5 <- NA
# dta.all$PROCHOICE5[dta.all$POS4_5=="Abortion should be available to anyone who wants it."] <- 4
# dta.all$PROCHOICE5[dta.all$POS4_5=="Abortion should be available, but with stricter limits than it is now."] <- 3
# dta.all$PROCHOICE5[dta.all$POS4_5=="Abortion should not be permitted except in cases of rape, incest, or when the life of the woman is at risk."] <- 2
# dta.all$PROCHOICE5[dta.all$POS4_5=="Abortion should not be permitted under any circumstances."] <- 1
# 
# dta.all$STAYIRAQ4 <- NA
# dta.all$STAYIRAQ4[dta.all$POS2_4=="The US should withdraw all troops from Iraq as soon as possible, regardless of conditions in Iraq."] <- 1
# dta.all$STAYIRAQ4[dta.all$POS2_4=="The US should set a deadline for withdrawing its troops if the Iraqi government doesn't show definite progress in traini"] <- 2
# dta.all$STAYIRAQ4[dta.all$POS2_4=="The US should keep its troops in Iraq as long as is needed until a stable democratic government is established there."] <- 3
# 
# ####
# 
# dta.all$TVNEWS6 <- NA
# dta.all$TVNEWS6[dta.all$ME1_1_7=="Yes" & ! dta.all$ME1_1_7 %in% c(NA)] <- 1
# dta.all$TVNEWS6[dta.all$ME1_1_7=="No" & ! dta.all$ME1_1_7 %in% c(NA)] <- 0
# 
# dta.all$NEWSPAPER6 <- NA
# dta.all$NEWSPAPER6[dta.all$ME1_2_7=="Yes" & ! dta.all$ME1_2_7 %in% c(NA)] <- 1
# dta.all$NEWSPAPER6[dta.all$ME1_2_7=="No" & ! dta.all$ME1_2_7 %in% c(NA)] <- 0
# 
# dta.all$TVTALK6 <- NA
# dta.all$TVTALK6[dta.all$ME1_3_7=="Yes" & ! dta.all$ME1_3_7 %in% c(NA)] <- 1
# dta.all$TVTALK6[dta.all$ME1_3_7=="No" & ! dta.all$ME1_3_7 %in% c(NA)] <- 0
# 
# dta.all$INTERNET6 <- NA
# dta.all$INTERNET6[dta.all$ME1_4_7=="Yes" & ! dta.all$ME1_4_7 %in% c(NA)] <- 1
# dta.all$INTERNET6[dta.all$ME1_4_7=="No" & ! dta.all$ME1_4_7 %in% c(NA)] <- 0
# 
# dta.all$RADIO6 <- NA
# dta.all$RADIO6[dta.all$ME1_5_7=="Yes" & ! dta.all$ME1_5_7 %in% c(NA)] <- 1
# dta.all$RADIO6[dta.all$ME1_5_7=="No" & ! dta.all$ME1_5_7 %in% c(NA)] <- 0
# 
# dta.all$MAGAZINE6 <- NA
# dta.all$MAGAZINE6[dta.all$ME1_6_7=="Yes" & ! dta.all$ME1_6_7 %in% c(NA)] <- 1
# dta.all$MAGAZINE6[dta.all$ME1_6_7=="No" & ! dta.all$ME1_6_7 %in% c(NA)] <- 0
# 
# dta.all$NEWSMEDIA6 <- dta.all$TVNEWS6+dta.all$NEWSPAPER6+dta.all$TVTALK6+dta.all$INTERNET6+dta.all$RADIO6+dta.all$MAGAZINE6
# 
# ##
# dta.all$TVNEWS10 <- NA
# dta.all$TVNEWS10[dta.all$ME1_1_10=="Yes" & ! dta.all$ME1_1_10 %in% c(NA)] <- 1
# dta.all$TVNEWS10[dta.all$ME1_1_10=="No" & ! dta.all$ME1_1_10 %in% c(NA)] <- 0
# 
# dta.all$NEWSPAPER10 <- NA
# dta.all$NEWSPAPER10[dta.all$ME1_2_10=="Yes" & ! dta.all$ME1_2_10 %in% c(NA)] <- 1
# dta.all$NEWSPAPER10[dta.all$ME1_2_10=="No" & ! dta.all$ME1_2_10 %in% c(NA)] <- 0
# 
# dta.all$TVTALK10 <- NA
# dta.all$TVTALK10[dta.all$ME1_3_10=="Yes" & ! dta.all$ME1_3_10 %in% c(NA)] <- 1
# dta.all$TVTALK10[dta.all$ME1_3_10=="No" & ! dta.all$ME1_3_10 %in% c(NA)] <- 0
# 
# dta.all$INTERNET10 <- NA
# dta.all$INTERNET10[dta.all$ME1_4_10=="Yes" & ! dta.all$ME1_4_10 %in% c(NA)] <- 1
# dta.all$INTERNET10[dta.all$ME1_4_10=="No" & ! dta.all$ME1_4_10 %in% c(NA)] <- 0
# 
# dta.all$RADIO10 <- NA
# dta.all$RADIO10[dta.all$ME1_5_10=="Yes" & ! dta.all$ME1_5_10 %in% c(NA)] <- 1
# dta.all$RADIO10[dta.all$ME1_5_10=="No" & ! dta.all$ME1_5_10 %in% c(NA)] <- 0
# 
# dta.all$MAGAZINE10 <- NA
# dta.all$MAGAZINE10[dta.all$ME1_6_10=="Yes" & ! dta.all$ME1_6_10 %in% c(NA)] <- 1
# dta.all$MAGAZINE10[dta.all$ME1_6_10=="No" & ! dta.all$ME1_6_10 %in% c(NA)] <- 0
# 
# dta.all$NEWSMEDIA10 <- dta.all$TVNEWS10+dta.all$NEWSPAPER10+dta.all$TVTALK10+dta.all$INTERNET10+dta.all$RADIO10+dta.all$MAGAZINE10
# cor(dta.all$NEWSMEDIA10,dta.all$NEWSMEDIA6,use="pairwise.complete.obs")
# #[1] 0.5629073
# 
# dta.all$EDYEARS4 <- NA
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="No formal education"] <- 0
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="1st, 2nd, 3rd, or 4th grade"] <- 4
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="5th or 6th grade"] <- 6
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="7th or 8th grade"] <- 8
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="9th grade"] <- 9
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="10th grade"] <- 10
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="11th grade"] <- 11
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="12th grade NO DIPLOMA"] <- 11.5
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="HIGH SCHOOL GRADUATE - high school DIPLOMA or the equivalent (GED)"] <- 12
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="Some college, no degree"] <- 13
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="Associate degree"] <- 14
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="Bachelor's degree"] <- 16
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="Master's degree"] <- 18
# dta.all$EDYEARS4[dta.all$PPEDUC_4=="Professional or Doctorate degree"] <- 19
# 
# dta.all$EDYEARS10 <- NA
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="No formal education"] <- 0
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="1st, 2nd, 3rd, or 4th grade"] <- 4
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="5th or 6th grade"] <- 6
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="7th or 8th grade"] <- 8
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="9th grade"] <- 9
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="10th grade"] <- 10
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="11th grade"] <- 11
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="12th grade NO DIPLOMA"] <- 11.5
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="HIGH SCHOOL GRADUATE - high school DIPLOMA or the equivalent (GED)"] <- 12
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="Some college, no degree"] <- 13
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="Associate degree"] <- 14
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="Bachelor's degree"] <- 16
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="Master's degree"] <- 18
# dta.all$EDYEARS10[dta.all$PPEDUC_10=="Professional or Doctorate degree"] <- 19
# 
# dta.all$TRUMPVC1012B <- apply(cbind(dta.all$TRUMPVC10B,dta.all$TRUMPVC11B,dta.all$TRUMPVC12B),1,mean)
# dta.all$ROMNEYVO67B <- apply(cbind(dta.all$ROMNEYVO6B,dta.all$ROMNEYVO7B),1,mean)
# dta.all$WBPREJUDICE1012 <- apply(cbind(dta.all$WBPREJUDICE10,dta.all$WBPREJUDICE11,dta.all$WBPREJUDICE12),1,mean)
# dta.all$WHPREJUDICE1012 <- apply(cbind(dta.all$WHPREJUDICE10,dta.all$WHPREJUDICE11,dta.all$WHPREJUDICE12),1,mean)
# 
# dta.all$HASHS4 <- 1*(dta.all$EDYEARS4 >= 12)
# dta.all$HASBA4 <- 1*(dta.all$EDYEARS4 >= 16)
# 
# dta.all$HASHS6 <- 1*(dta.all$EDYEARS6 >= 12)
# dta.all$HASBA6 <- 1*(dta.all$EDYEARS6 >= 16)
# dta.all$OVER656 <- 1*(dta.all$AGE6 >= 65)
# 
# dta.all$HASHS10 <- 1*(dta.all$EDYEARS10 >= 12)
# dta.all$HASBA10 <- 1*(dta.all$EDYEARS10 >= 16)
# 
# dta.all$AGE4 <- dta.all$PPAGE_4
# dta.all$OVER65 <- 1*(dta.all$PPAGE_4 >= 65)
# 
# dta.all$HISP4 <- 1*(dta.all$PPETHM_4=="Hispanic")
# dta.all$BLACK4 <- 1*(dta.all$PPETHM_4=="Black, Non-Hispanic")
# dta.all$WHITE4 <- 1*(dta.all$PPETHM_4=="White, Non-Hispanic")
# 
# dta.all$MCCAINVO5B <- NA
# dta.all$MCCAINVO5B[dta.all$MCCAINVO5A==0] <- -1
# dta.all$MCCAINVO5B[dta.all$MCCAINVO5A==9] <- 0
# dta.all$MCCAINVO5B[dta.all$MCCAINVO5A==1] <- 1
# 
# dta.all$RUBIOVC10B <- NA
# dta.all$RUBIOVC10B[dta.all$RUBIOVC10A==1] <- 1
# dta.all$RUBIOVC10B[dta.all$RUBIOVC10A==0] <- -1
# dta.all$RUBIOVC10B[dta.all$RUBIOVC10A==9] <- 0

# dta.all$CRUZVC10B <- NA
# dta.all$CRUZVC10B[dta.all$CRUZVC10A==1] <- 1
# dta.all$CRUZVC10B[dta.all$CRUZVC10A==0] <- -1
# dta.all$CRUZVC10B[dta.all$CRUZVC10A==9] <- 0
# 
# dta.all$CLINTON10P <- dta.all$CLINTOND10

### subset data for public release----
vars1 <- unique(c("INCOME1","EDYEARS4","HASHS4","HASBA4","PID4","PID0","UNION0","CATHOLIC0","PROTESTANT0","FEMALE6",
                  "AGE4","OVER65","BLACK4","HISP4","WHITE4","BLACK6","HISP6","WHITE6",
                  "WBPREJUDICE4R","WBPREJUDICE5R","WBPREJUDICE6R","WBPREJUDICE7R","WBPREJUDICE10R","WBPREJUDICE11R","WBPREJUDICE12R",
                  "WHPREJUDICE6R","WHPREJUDICE7R","WHPREJUDICE10R","WHPREJUDICE11R","WBPREJUDICE67","WBPREJUDICE45",
                  
                  "GOPVDEM06B","WBPREJUDICE45","INCOME1F","UNION0","CATHOLIC0","PROTESTANT0","EDYEARS6","FEMALE6","AGE6","BUSHVK0A1",
                  "BUSHVK0A9",
                  "WEAKDEM0","INDDEM0","INDEP0","INDGOP0","WEAKGOP0","STRONGGOP0",
                  "WEAKDEM6","INDDEM6","INDEP6","INDGOP6","WEAKGOP6","STRONGGOP6","PID6",
                  
                  "MCCAINVO3B","MCCAINVO4A1","MCCAINVO4A9","MCCAINVO4B","CRUZVC10B","ROMNEYVO6B","ROMNEYVO6A1","ROMNEYVO6A9","ROMNEYVO7B","RUBIOVC10B",
                  "MCCAINVO5B","TRUMPVC10B","TRUMPVC11B","TRUMPVC12B","TRUMPVC10A","TRUMPVC12A","GOPREPVD14B","TRUMPVC11F","TRUMPVC12F",
                  "WBPREJUDICE4", "WBPREJUDICE67", "WHPREJUDICE67", "INCOME1F","UNION0" ,"CATHOLIC0",
                  "PROTESTANT0",   "EDYEARS6" , "FEMALE6","AGE6",
                  "TRUMP10P","CLINTON10P","SANDERS10P","VP1A_REP_1","VP1A_DEM_1"))

dta.sub <- subset(dta.all,select=vars1)
save(dta.sub,file="voting-pb-replication-hopkins-07312019.Rdata")

### start here for replication: load & subset data ----
load("voting-pb-replication-hopkins-07312019.Rdata")

dta.finalwave <- dta.sub[! dta.sub$TRUMPVC12A %in% c(NA),]
dta.nhw.final <- dta.nhw <- dta.sub[dta.sub$WHITE6==1 & ! dta.sub$TRUMPVC12A %in% c(NA),]
rownames(dta.nhw) <- 1:dim(dta.nhw)[1]
dta.nhw2 <- dta.sub[dta.sub$WHITE6==1 & ! dta.sub$TRUMPVC10A %in% c(NA),]
rownames(dta.nhw2) <- 1:dim(dta.nhw2)[1]

#### DESCRIPTIVE STATISTICS -----

###appendix table 2; partial coverage
vars <- c("INCOME1","EDYEARS4","HASHS4","HASBA4","PID4","UNION0","CATHOLIC0","PROTESTANT0","FEMALE6",
          "AGE4","OVER65","BLACK4","HISP4","WHITE4")
rmat <- matrix(NA,length(vars),5)
for(i in 1:length(vars)){
  txt <- paste("hold <- dta.finalwave$",vars[i])
  eval(parse(text=txt))
  rmat[i,1] <- min(hold,na.rm=T)
  rmat[i,2] <- mean(hold,na.rm=T)
  rmat[i,3] <- max(hold,na.rm=T)
  rmat[i,4] <- sum(hold %in% c(NA))/length(hold)
}
rownames(rmat) <- vars
colnames(rmat) <- c("Min","Mean","Max","Pct Missing","Census")
xtable(rmat,digits=c(0,2,2,2,2,2))

### non-hispanic white descriptives
### Appendix Table 3
vars <- c("INCOME1","EDYEARS4","HASHS4","HASBA4","PID4","UNION0","CATHOLIC0","PROTESTANT0","FEMALE6","AGE4","OVER65","BLACK6","HISP6","WHITE6")
rmat <- matrix(NA,length(vars),4)
for(i in 1:length(vars)){
  
  txt <- paste("hold <- dta.nhw$",vars[i])
  eval(parse(text=txt))
  
  rmat[i,1] <- min(hold,na.rm=T)
  rmat[i,2] <- mean(hold,na.rm=T)
  rmat[i,3] <- max(hold,na.rm=T)
  rmat[i,4] <- sum(hold %in% c(NA))/length(hold)
  
}
rownames(rmat) <- vars
colnames(rmat) <- c("Min","Mean","Max","Pct Missing")
xtable(rmat,digits=c(0,2,2,2,2))

### Appendix Table 4
vars <- c("WBPREJUDICE4R","WBPREJUDICE5R","WBPREJUDICE6R","WBPREJUDICE7R","WBPREJUDICE10R","WBPREJUDICE11R","WBPREJUDICE12R")
rmat <- matrix(NA,length(vars),5)
for(i in 1:length(vars)){
	
	txt <- paste("hold <- dta.nhw$",vars[i])
	eval(parse(text=txt))

	rmat[i,1] <- min(hold,na.rm=T)
	rmat[i,2] <- mean(hold,na.rm=T)
	rmat[i,3] <- max(hold,na.rm=T)
	rmat[i,4] <- sd(hold,na.rm=T)
	rmat[i,5] <- sum(hold %in% c(NA))/length(hold)
	
}
rownames(rmat) <- vars
colnames(rmat) <- c("Min","Mean","Max","SD","Pct Missing")
xtable(rmat,digits=c(0,3,3,3,3,3))

### Appendix Table 5
vars <- c("WHPREJUDICE6R","WHPREJUDICE7R","WHPREJUDICE10R","WHPREJUDICE11R")
rmat <- matrix(NA,length(vars),5)
for(i in 1:length(vars)){
	
	txt <- paste("hold <- dta.nhw.final$",vars[i])
	eval(parse(text=txt))

	rmat[i,1] <- min(hold,na.rm=T)
	rmat[i,2] <- mean(hold,na.rm=T)
	rmat[i,3] <- max(hold,na.rm=T)
	rmat[i,4] <- sd(hold,na.rm=T)
	rmat[i,5] <- sum(hold %in% c(NA))/length(hold)
	
}
rownames(rmat) <- vars
colnames(rmat) <- c("Min","Mean","Max","SD","Pct Missing")
xtable(rmat,digits=c(0,3,3,3,3,3))

### MNP estimates for anti-Black prejudice -----
res0c <- mnp(GOPVDEM06B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw,n.draws=250)
lout0c <- lm(GOPVDEM06B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw)

mf <- model.frame(lout0c)
colnames(mf)[1] <- "GOPVDEM06B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.2,na.rm=T)
mf2$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.8,na.rm=T)

pout0cb <- predict.mnp(res0c,newdata=mf1)
zero.0cb <- apply(apply(pout0cb$y,2,is.zero),2,mean)
quantile(zero.0cb,c(0.025,.5,.975))

one.0cb <- apply(apply(pout0cb$y,2,is.one),2,mean)
quantile(one.0cb,c(0.025,.5,.975))

two.0cb <- apply(apply(pout0cb$y,2,is.two),2,mean)
quantile(two.0cb,c(0.025,.5,.975))

pout0cc <- predict.mnp(res0c,newdata=mf2)
zero.0cc <- apply(apply(pout0cc$y,2,is.zero),2,mean)
quantile(zero.0c,c(0.025,.5,.975))

one.0cc <- apply(apply(pout0cc$y,2,is.one),2,mean)
quantile(one.0cc,c(0.025,.5,.975))

two.0cc <- apply(apply(pout0cc$y,2,is.two),2,mean)
quantile(two.0cc,c(0.025,.5,.975))

gop.diff.0cc <- two.0cc-two.0cb
quantile(gop.diff.0cc,c(0.025,.5,.975))

###

res3 <- mnp(MCCAINVO3B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw,n.draws=5000)
lout3 <- lm(MCCAINVO3B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw)

mf <- model.frame(lout3)
colnames(mf)[1] <- "MCCAINVO3B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.2,na.rm=T)
mf2$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.8,na.rm=T)

pout3b <- predict.mnp(res3,newdata=mf1)
zero.3b <- apply(apply(pout3b$y,2,is.zero),2,mean)
quantile(zero.3b,c(0.025,.5,.975))

one.3b <- apply(apply(pout3b$y,2,is.one),2,mean)
quantile(one.3b,c(0.025,.5,.975))

two.3b <- apply(apply(pout3b$y,2,is.two),2,mean)
quantile(two.3b,c(0.025,.5,.975))

pout3c <- predict.mnp(res3,newdata=mf2)

zero.3c <- apply(apply(pout3c$y,2,is.zero),2,mean)
quantile(zero.3c,c(0.025,.5,.975))

one.3c <- apply(apply(pout3c$y,2,is.one),2,mean)
quantile(one.3c,c(0.025,.5,.975))

two.3c <- apply(apply(pout3c$y,2,is.two),2,mean)
quantile(two.3c,c(0.025,.5,.975))

gop.diff.3 <- two.3c-two.3b
quantile(gop.diff.3,c(0.025,.5,.975))

###
res4 <- mnp(MCCAINVO4B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw,n.draws=4000)
lout4 <- lm(MCCAINVO4B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw)

mf <- model.frame(lout4)
colnames(mf)[1] <- "MCCAINVO4B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.2,na.rm=T)
mf2$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.8,na.rm=T)

pout4b <- predict.mnp(res4,newdata=mf1)
zero.4b <- apply(apply(pout4b$y,2,is.zero),2,mean)
quantile(zero.4b,c(0.025,.5,.975))

one.4b <- apply(apply(pout4b$y,2,is.one),2,mean)
quantile(one.4b,c(0.025,.5,.975))

two.4b <- apply(apply(pout4b$y,2,is.two),2,mean)
quantile(two.4b,c(0.025,.5,.975))

pout4c <- predict.mnp(res4,newdata=mf2)
zero.4c <- apply(apply(pout4c$y,2,is.zero),2,mean)
quantile(zero.4c,c(0.025,.5,.975))

one.4c <- apply(apply(pout4c$y,2,is.one),2,mean)
quantile(one.4c,c(0.025,.5,.975))

two.4c <- apply(apply(pout4c$y,2,is.two),2,mean)
quantile(two.4c,c(0.025,.5,.975))

gop.diff.4 <- two.4c-two.4b
quantile(gop.diff.4,c(0.025,.5,.975))


res5 <- mnp(MCCAINVO5B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw,n.draws=2500)
lout5 <- lm(MCCAINVO5B ~ WBPREJUDICE45+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+BUSHVK0A1+BUSHVK0A9+WEAKDEM0+INDDEM0+INDEP0+INDGOP0+WEAKGOP0+STRONGGOP0,data=dta.nhw)

mf <- model.frame(lout5)
colnames(mf)[1] <- "MCCAINVO5B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.2,na.rm=T)
mf2$WBPREJUDICE45 <- quantile(dta.nhw$WBPREJUDICE45,.8,na.rm=T)

pout5b <- predict.mnp(res5,newdata=mf1)
zero.5b <- apply(apply(pout5b$y,2,is.zero),2,mean)
quantile(zero.5b,c(0.025,.5,.975))

one.5b <- apply(apply(pout5b$y,2,is.one),2,mean)
quantile(one.5b,c(0.025,.5,.975))

two.5b <- apply(apply(pout5b$y,2,is.two),2,mean)
quantile(two.5b,c(0.025,.5,.975))

pout5c <- predict.mnp(res5,newdata=mf2)

zero.5c <- apply(apply(pout5c$y,2,is.zero),2,mean)
quantile(zero.5c,c(0.025,.5,.975))

one.5c <- apply(apply(pout5c$y,2,is.one),2,mean)
quantile(one.5c,c(0.025,.5,.975))

two.5c <- apply(apply(pout5c$y,2,is.two),2,mean)
quantile(two.5c,c(0.025,.5,.975))

gop.diff.5 <- two.5c-two.5b
quantile(gop.diff.5,c(0.025,.5,.975))

res6 <- mnp(ROMNEYVO6B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout6 <- lm(ROMNEYVO6B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res6)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res6)$coef.table[s1,]
t2 <- summary(res6)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Intercept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
xtable(tjoint,digits=c(0,3,3,3,3))

mf <- model.frame(lout6)
colnames(mf)[1] <- "ROMNEYVO6B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout6b <- predict.mnp(res6,newdata=mf1)
zero.6b <- apply(apply(pout6b$y,2,is.zero),2,mean)
quantile(zero.6b,c(0.025,.5,.975))

one.6b <- apply(apply(pout6b$y,2,is.one),2,mean)
quantile(one.6b,c(0.025,.5,.975))

two.6b <- apply(apply(pout6b$y,2,is.two),2,mean)
quantile(two.6b,c(0.025,.5,.975))

#table(pout8b$y)/sum(table(pout8b$y))

pout6c <- predict.mnp(res6,newdata=mf2)
zero.6c <- apply(apply(pout6c$y,2,is.zero),2,mean)
quantile(zero.6c,c(0.025,.5,.975))

one.6c <- apply(apply(pout6c$y,2,is.one),2,mean)
quantile(one.6c,c(0.025,.5,.975))

two.6c <- apply(apply(pout6c$y,2,is.two),2,mean)
quantile(two.6c,c(0.025,.5,.975))

gop.diff.6 <- two.6c-two.6b
quantile(gop.diff.6,c(0.025,.5,.975))

####
res7 <- mnp(ROMNEYVO7B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=2000)#,trace=T)#,invcdf=T)
lout7 <- lm(ROMNEYVO7B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res7)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res7)$coef.table[s1,]
t2 <- summary(res7)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Intercept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
xtable(tjoint,digits=c(0,3,3,3,3))

mf <- model.frame(lout7)
colnames(mf)[1] <- "ROMNEYVO7B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout7b <- predict.mnp(res7,newdata=mf1)
zero.7b <- apply(apply(pout7b$y,2,is.zero),2,mean)
quantile(zero.7b,c(0.025,.5,.975))

one.7b <- apply(apply(pout6b$y,2,is.one),2,mean)
quantile(one.7b,c(0.025,.5,.975))

two.7b <- apply(apply(pout7b$y,2,is.two),2,mean)
quantile(two.7b,c(0.025,.5,.975))

pout7c <- predict.mnp(res7,newdata=mf2)
zero.7c <- apply(apply(pout7c$y,2,is.zero),2,mean)
quantile(zero.7c,c(0.025,.5,.975))

one.7c <- apply(apply(pout7c$y,2,is.one),2,mean)
quantile(one.7c,c(0.025,.5,.975))

two.7c <- apply(apply(pout7c$y,2,is.two),2,mean)
quantile(two.7c,c(0.025,.5,.975))

gop.diff.7 <- two.7c-two.7b
quantile(gop.diff.7,c(0.025,.5,.975))

res8 <- mnp(GOPREPVD14B~WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout8 <- lm(GOPREPVD14B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res8)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res8)$coef.table[s1,]
t2 <- summary(res8)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Intercept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
xtable(tjoint,digits=c(0,3,3,3,3))


mf <- model.frame(lout8)
colnames(mf)[1] <- "GOPREPVD14B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout8b <- predict.mnp(res8,newdata=mf1)
zero.8b <- apply(apply(pout8b$y,2,is.zero),2,mean)
quantile(zero.8b,c(0.025,.5,.975))

one.8b <- apply(apply(pout8b$y,2,is.one),2,mean)
quantile(one.8b,c(0.025,.5,.975))

two.8b <- apply(apply(pout8b$y,2,is.two),2,mean)
quantile(two.8b,c(0.025,.5,.975))

pout8c <- predict.mnp(res8,newdata=mf2)
zero.8c <- apply(apply(pout8c$y,2,is.zero),2,mean)
quantile(zero.8c,c(0.025,.5,.975))

one.8c <- apply(apply(pout8c$y,2,is.one),2,mean)
quantile(one.8c,c(0.025,.5,.975))

two.8c <- apply(apply(pout8c$y,2,is.two),2,mean)
quantile(two.8c,c(0.025,.5,.975))

gop.diff.8 <- two.8c-two.8b
quantile(gop.diff.8,c(0.025,.5,.975))

res10 <- mnp(TRUMPVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=2500)
lout10 <- lm(TRUMPVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res10)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res10)$coef.table[s1,]
t2 <- summary(res10)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Intercept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
### Appendix Table 10 (may not match due to sampling)
xtable(tjoint,digits=c(0,3,3,3,3))

mf <- model.frame(lout10)
colnames(mf)[1] <- "TRUMPVC10B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout10b <- predict.mnp(res10,newdata=mf1)
zero.10b <- apply(apply(pout10b$y,2,is.zero),2,mean)
quantile(zero.10b,c(0.025,.5,.975))

one.10b <- apply(apply(pout10b$y,2,is.one),2,mean)
quantile(one.10b,c(0.025,.5,.975))

two.10b <- apply(apply(pout10b$y,2,is.two),2,mean)
quantile(two.10b,c(0.025,.5,.975))

pout10c <- predict.mnp(res10,newdata=mf2)
zero.10c <- apply(apply(pout10c$y,2,is.zero),2,mean)
quantile(zero.10c,c(0.025,.5,.975))

one.10c <- apply(apply(pout10c$y,2,is.one),2,mean)
quantile(one.10c,c(0.025,.5,.975))

two.10c <- apply(apply(pout10c$y,2,is.two),2,mean)
quantile(two.10c,c(0.025,.5,.975))

gop.diff.10.bp <- gop.diff.10 <- two.10c-two.10b
quantile(gop.diff.10,c(0.025,.5,.975))
1-sum(gop.diff.10 > 0)/length(gop.diff.10)

### 
gop.diff.10.6 <- gop.diff.10-gop.diff.6
1-sum(gop.diff.10.6 >0)/length(gop.diff.10.6)

####
res11 <- mnp(TRUMPVC11B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout11 <- lm(TRUMPVC11B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res11)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res11)$coef.table[s1,]
t2 <- summary(res11)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Interecept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
### Appendix Table 11 (may not exactly match due to simulation differences)
xtable(tjoint,digits=c(0,3,3,3,3))

mf <- model.frame(lout11)
#mf <- model.frame(lout10)
colnames(mf)[1] <- "TRUMPVC11B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout11b <- predict.mnp(res11,newdata=mf1)
zero.11b <- apply(apply(pout11b$y,2,is.zero),2,mean)
quantile(zero.11b,c(0.025,.5,.975))

one.11b <- apply(apply(pout10b$y,2,is.one),2,mean)
quantile(one.11b,c(0.025,.5,.975))

two.11b <- apply(apply(pout11b$y,2,is.two),2,mean)
quantile(two.11b,c(0.025,.5,.975))

pout11c <- predict.mnp(res11,newdata=mf2)
zero.11c <- apply(apply(pout11c$y,2,is.zero),2,mean)
quantile(zero.11c,c(0.025,.5,.975))

one.11c <- apply(apply(pout11c$y,2,is.one),2,mean)
quantile(one.11c,c(0.025,.5,.975))

two.11c <- apply(apply(pout11c$y,2,is.two),2,mean)
quantile(two.11c,c(0.025,.5,.975))

gop.diff.11 <- two.11c-two.11b
quantile(gop.diff.11,c(0.025,.5,.975))

gop.diff.11.6 <- gop.diff.11-gop.diff.6
sum(gop.diff.11.6 >0)/length(gop.diff.11.6)

###
res12 <- mnp(TRUMPVC12F ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout12 <- lm(TRUMPVC12B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

n.param <- dim(summary(res12)$coef.table)[1]
s1 <- seq(from=1,to=(n.param)-1,by=2)
s2 <- seq(from=2,to=(n.param),by=2)

t1 <- summary(res12)$coef.table[s1,]
t2 <- summary(res12)$coef.table[s2,]

tjoint <- cbind(t1[,1:2],t2[,1:2])

rn.new <- c("Intercept","Anti-Black Prejudice '12","Anti-Latino Prejudice '12", 
	"Income '08","Union '08","Catholic '08","Protestant '08","Education '08",
	"Female '08","Age '08","Lagged GOP Support","Lagged Support Neither","Weak Dem. '08",
	"Lean Dem. '08","Independent '08","Lean GOP '08","Weak GOP '08","Strong GOP '08")

cbind(rn.new,rownames(tjoint))
rownames(tjoint) <- rn.new
### Appendix Table 12; may not match exactly due to rounding
xtable(tjoint,digits=c(0,3,3,3,3))

mf <- model.frame(lout12)
colnames(mf)[1] <- "TRUMPVC12F"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout12b <- predict.mnp(res12,newdata=mf1)
zero.12b <- apply(apply(pout12b$y,2,is.zero),2,mean)
quantile(zero.12b,c(0.025,.5,.975))

one.12b <- apply(apply(pout12b$y,2,is.one),2,mean)
quantile(one.12b,c(0.025,.5,.975))

two.12b <- apply(apply(pout12b$y,2,is.two),2,mean)
quantile(two.12b,c(0.025,.5,.975))

pout12c <- predict.mnp(res12,newdata=mf2)
zero.12c <- apply(apply(pout12c$y,2,is.zero),2,mean)
quantile(zero.12c,c(0.025,.5,.975))

one.12c <- apply(apply(pout12c$y,2,is.one),2,mean)
quantile(one.12c,c(0.025,.5,.975))

two.12c <- apply(apply(pout12c$y,2,is.two),2,mean)
quantile(two.12c,c(0.025,.5,.975))

gop.diff.12 <- two.12c-two.12b
quantile(gop.diff.12,c(0.025,.5,.975))

gop.diff.12.6 <- gop.diff.12-gop.diff.6
sum(gop.diff.12.6 >0)/length(gop.diff.12.6)
#[1] 0.8301

rmat <- as.data.frame(matrix(NA,8,2))

rmat[1,1] <- mean(gop.diff.0cc)
rmat[1,2] <- sd(gop.diff.0cc)

rmat[2,1] <- mean(gop.diff.4)
rmat[2,2] <- sd(gop.diff.4)

rmat[3,1] <- mean(gop.diff.5)
rmat[3,2] <- sd(gop.diff.5)

rmat[4,1] <- mean(gop.diff.6)
rmat[4,2] <- sd(gop.diff.6)

rmat[5,1] <- mean(gop.diff.7)
rmat[5,2] <- sd(gop.diff.7)

rmat[6,1] <- mean(gop.diff.8)
rmat[6,2] <- sd(gop.diff.8)

rmat[7,1] <- mean(gop.diff.10)
rmat[7,2] <- sd(gop.diff.10)

rmat[8,1] <- mean(gop.diff.11)
rmat[8,2] <- sd(gop.diff.11)

rmat[9,1] <- mean(gop.diff.12)
rmat[9,2] <- sd(gop.diff.12)

rmat$date <- c(2006.9166667,2008.8333,2008.9166667,2012.8333,2013,2014.8333333,2016.083333,2016.83333,2017)

xtable(rmat[-c(1:3),c(1:2)],digits=c(0,3,3))

vars <- c("gop.diff.6","gop.diff.7","gop.diff.8","gop.diff.10","gop.diff.11","gop.diff.12")
pmat <- matrix(NA,length(vars),length(vars))
for(i in 1:length(vars)){
	for(j in 1:length(vars)){
		txt1 <- paste("var1 <- ",vars[i])
		eval(parse(text=txt1))
		txt2 <- paste("var2 <- ",vars[j])
		eval(parse(text=txt2))
		n.min <- min(c(length(var1),length(var2)))
		pmat[i,j] <- sum(var1[1:n.min] > var2[1:n.min])/n.min
	}
}
xtable(pmat,digits=c(0,rep(2,length(vars))))

rmat$id <- c(rep(1,6),rep(2,3))
rmat$var1 <- rmat[,1]
rmat$SE1 <- rmat[,2]

library(plyr)
#f1()
ddply(rmat, "id", summarize, 
      newMean = weighted.mean(x=var1, 1/SE1, na.rm = TRUE),
      newSE = 1/sum(1/SE1, na.rm = TRUE)
      )
      
#  id    newMean       newSE
#1  1 0.01288014 0.004221450
#2  2 0.04752158 0.007879439      
      
##### figure 1
pdf("lenz-test-mnp-noarrows-2012-2016-02112018.pdf",width=10)

plot(x=rmat$date[-c(1:3)],y=rmat[-c(1:3),1],pch=16,ylim=c(-.1,.10),xlim=c(2012,2017),cex.main=2,main="Change in GOP Support \n as Anti-Black Prejudice Rises",ylab="Shift in Prob. GOP Support",xlab="Date",cex.lab=1.3,xaxt="n")
axis(side=1,at=c(2012,2013,2014,2015,2016,2017),labels=c("Jan12","Jan13","Jan14","Jan15","Jan16","Jan17"))

for(i in 4:dim(rmat)[1]){
	lines(x=c(rmat$date[i],rmat$date[i]),y=c(rmat[i,1]-rmat[i,2]*1.96,rmat[i,1]+rmat[i,2]*1.96))	
}
text(x=2012.8333,y=-.09,"Romney \n vs. Obama")
text(x=2014.8333333,y=-.09,"House GOP \n vs. Dem.")
text(x=2016.5,y=-.09,"Trump \n vs. Clinton")

abline(h=0)

dev.off()






### MNP Estimates for Latino prejudice ----

res6 <- mnp(ROMNEYVO6B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=5000)
lout6 <- lm(ROMNEYVO6B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout6)
colnames(mf)[1] <- "ROMNEYVO6B"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout6b <- predict.mnp(res6,newdata=mf1)
zero.6b <- apply(apply(pout6b$y,2,is.zero),2,mean)
quantile(zero.6b,c(0.025,.5,.975))

one.6b <- apply(apply(pout6b$y,2,is.one),2,mean)
quantile(one.6b,c(0.025,.5,.975))

two.6b <- apply(apply(pout6b$y,2,is.two),2,mean)
quantile(two.6b,c(0.025,.5,.975))

pout6c <- predict.mnp(res6,newdata=mf2)
zero.6c <- apply(apply(pout6c$y,2,is.zero),2,mean)
quantile(zero.6c,c(0.025,.5,.975))

one.6c <- apply(apply(pout6c$y,2,is.one),2,mean)
quantile(one.6c,c(0.025,.5,.975))

two.6c <- apply(apply(pout6c$y,2,is.two),2,mean)
quantile(two.6c,c(0.025,.5,.975))

gop.diff.6 <- two.6c-two.6b
quantile(gop.diff.6,c(0.025,.5,.975))

res7 <- mnp(ROMNEYVO7B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=2000)#,trace=T)#,invcdf=T)
lout7 <- lm(ROMNEYVO7B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+MCCAINVO4A1+MCCAINVO4A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout7)
colnames(mf)[1] <- "ROMNEYVO7B"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout7b <- predict.mnp(res7,newdata=mf1)
zero.7b <- apply(apply(pout7b$y,2,is.zero),2,mean)
quantile(zero.7b,c(0.025,.5,.975))

one.7b <- apply(apply(pout6b$y,2,is.one),2,mean)
quantile(one.7b,c(0.025,.5,.975))

two.7b <- apply(apply(pout7b$y,2,is.two),2,mean)
quantile(two.7b,c(0.025,.5,.975))

pout7c <- predict.mnp(res7,newdata=mf2)
zero.7c <- apply(apply(pout7c$y,2,is.zero),2,mean)
quantile(zero.7c,c(0.025,.5,.975))

one.7c <- apply(apply(pout7c$y,2,is.one),2,mean)
quantile(one.7c,c(0.025,.5,.975))

two.7c <- apply(apply(pout7c$y,2,is.two),2,mean)
quantile(two.7c,c(0.025,.5,.975))

gop.diff.7 <- two.7c-two.7b
quantile(gop.diff.7,c(0.025,.5,.975))

res8 <- mnp(GOPREPVD14B~WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout8 <- lm(GOPREPVD14B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout8)
colnames(mf)[1] <- "GOPREPVD14B"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout8b <- predict.mnp(res8,newdata=mf1)
zero.8b <- apply(apply(pout8b$y,2,is.zero),2,mean)
quantile(zero.8b,c(0.025,.5,.975))

one.8b <- apply(apply(pout8b$y,2,is.one),2,mean)
quantile(one.8b,c(0.025,.5,.975))

two.8b <- apply(apply(pout8b$y,2,is.two),2,mean)
quantile(two.8b,c(0.025,.5,.975))

pout8c <- predict.mnp(res8,newdata=mf2)
zero.8c <- apply(apply(pout8c$y,2,is.zero),2,mean)
quantile(zero.8c,c(0.025,.5,.975))

one.8c <- apply(apply(pout8c$y,2,is.one),2,mean)
quantile(one.8c,c(0.025,.5,.975))

two.8c <- apply(apply(pout8c$y,2,is.two),2,mean)
quantile(two.8c,c(0.025,.5,.975))

gop.diff.8 <- two.8c-two.8b
quantile(gop.diff.8,c(0.025,.5,.975))

res10 <- mnp(TRUMPVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=5000)
lout10 <- lm(TRUMPVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout10)
colnames(mf)[1] <- "TRUMPVC10B"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout10b <- predict.mnp(res10,newdata=mf1)
zero.10b <- apply(apply(pout10b$y,2,is.zero),2,mean)
quantile(zero.10b,c(0.025,.5,.975))

one.10b <- apply(apply(pout10b$y,2,is.one),2,mean)
quantile(one.10b,c(0.025,.5,.975))

two.10b <- apply(apply(pout10b$y,2,is.two),2,mean)
quantile(two.10b,c(0.025,.5,.975))

pout10c <- predict.mnp(res10,newdata=mf2)
zero.10c <- apply(apply(pout10c$y,2,is.zero),2,mean)
quantile(zero.10c,c(0.025,.5,.975))

one.10c <- apply(apply(pout10c$y,2,is.one),2,mean)
quantile(one.10c,c(0.025,.5,.975))

two.10c <- apply(apply(pout10c$y,2,is.two),2,mean)
quantile(two.10c,c(0.025,.5,.975))

gop.diff.10 <- two.10c-two.10b
quantile(gop.diff.10,c(0.025,.5,.975))
1-sum(gop.diff.10 > 0)/length(gop.diff.10)

gop.diff.10.6 <- gop.diff.10-gop.diff.6
1-sum(gop.diff.10.6 >0)/length(gop.diff.10.6)

res11 <- mnp(TRUMPVC11F ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout11 <- lm(TRUMPVC11B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout11)
colnames(mf)[1] <- "TRUMPVC11F"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout11b <- predict.mnp(res11,newdata=mf1)
zero.11b <- apply(apply(pout11b$y,2,is.zero),2,mean)
quantile(zero.11b,c(0.025,.5,.975))

one.11b <- apply(apply(pout10b$y,2,is.one),2,mean)
quantile(one.11b,c(0.025,.5,.975))

two.11b <- apply(apply(pout11b$y,2,is.two),2,mean)
quantile(two.11b,c(0.025,.5,.975))

pout11c <- predict.mnp(res11,newdata=mf2)
zero.11c <- apply(apply(pout11c$y,2,is.zero),2,mean)
quantile(zero.11c,c(0.025,.5,.975))

one.11c <- apply(apply(pout11c$y,2,is.one),2,mean)
quantile(one.11c,c(0.025,.5,.975))

two.11c <- apply(apply(pout11c$y,2,is.two),2,mean)
quantile(two.11c,c(0.025,.5,.975))

gop.diff.11 <- two.11c-two.11b
quantile(gop.diff.11,c(0.025,.5,.975))

gop.diff.11.6 <- gop.diff.11-gop.diff.6
sum(gop.diff.11.6 >0)/length(gop.diff.11.6)

res12 <- mnp(TRUMPVC12F ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=10000)
lout12 <- lm(TRUMPVC12B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout12)
colnames(mf)[1] <- "TRUMPVC12F"
mf1 <- mf2 <- mf
mf1$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.2,na.rm=T)
mf2$WHPREJUDICE67 <- quantile(dta.nhw$WHPREJUDICE67,.8,na.rm=T)

pout12b <- predict.mnp(res12,newdata=mf1)
zero.12b <- apply(apply(pout12b$y,2,is.zero),2,mean)
quantile(zero.12b,c(0.025,.5,.975))

one.12b <- apply(apply(pout12b$y,2,is.one),2,mean)
quantile(one.12b,c(0.025,.5,.975))

two.12b <- apply(apply(pout12b$y,2,is.two),2,mean)
quantile(two.12b,c(0.025,.5,.975))

pout12c <- predict.mnp(res12,newdata=mf2)
zero.12c <- apply(apply(pout12c$y,2,is.zero),2,mean)
quantile(zero.12c,c(0.025,.5,.975))

one.12c <- apply(apply(pout12c$y,2,is.one),2,mean)
quantile(one.12c,c(0.025,.5,.975))

two.12c <- apply(apply(pout12c$y,2,is.two),2,mean)
quantile(two.12c,c(0.025,.5,.975))

gop.diff.12 <- two.12c-two.12b
quantile(gop.diff.12,c(0.025,.5,.975))

gop.diff.12.6 <- gop.diff.12-gop.diff.6
sum(gop.diff.12.6 >0)/length(gop.diff.12.6)

rmath <- as.data.frame(matrix(NA,8,2))

rmath[4,1] <- mean(gop.diff.6)
rmath[4,2] <- sd(gop.diff.6)

rmath[5,1] <- mean(gop.diff.7)
rmath[5,2] <- sd(gop.diff.7)

rmath[6,1] <- mean(gop.diff.8)
rmath[6,2] <- sd(gop.diff.8)

rmath[7,1] <- mean(gop.diff.10)
rmath[7,2] <- sd(gop.diff.10)

rmath[8,1] <- mean(gop.diff.11)
rmath[8,2] <- sd(gop.diff.11)

rmath[9,1] <- mean(gop.diff.12)
rmath[9,2] <- sd(gop.diff.12)

rmath$date <- c(2006.9166667,2008.8333,2008.9166667,2012.8333,2013,2014.8333333,2016.083333,2016.83333,2017)

xtable(rmath[,c(1:2)],digits=c(0,3,3))

### make bottom panel of Figure 1
#pdf("lenz-test-hispanic-mnp-noarrows-2012-2016-02112018.pdf",width=10)
plot(x=rmath$date[-c(1:3)],y=rmath[-c(1:3),1],pch=16,ylim=c(-.1,.10),xlim=c(2012,2017),cex.main=2,main="Change in GOP Support \n as Anti-Latino Prejudice Rises",ylab="Shift in Prob. GOP Support",xlab="Date",cex.lab=1.3,xaxt="n")
axis(side=1,at=c(2012,2013,2014,2015,2016,2017),labels=c("Jan12","Jan13","Jan14","Jan15","Jan16","Jan17"))

for(i in 4:dim(rmath)[1]){
  lines(x=c(rmath$date[i],rmath$date[i]),y=c(rmath[i,1]-rmath[i,2]*1.96,rmath[i,1]+rmath[i,2]*1.96))	
}
text(x=2012.8333,y=-.09,"Romney \n vs. Obama")
text(x=2014.8333333,y=-.09,"House GOP \n vs. Dem.")
text(x=2016.5,y=-.09,"Trump \n vs. Clinton")
abline(h=0)
#dev.off()


### fit JOINT MODEL ----
mf6 <- model.frame(lout6)
mf7 <- model.frame(lout7)
mf10 <- model.frame(lout10)
mf11 <- model.frame(lout11)
mf12 <- model.frame(lout12)

colnames(mf12) <- colnames(mf11) <- colnames(mf10) <- colnames(mf7) <- colnames(mf6) <- c( "GOPVOTING",    "WBPREJUDICE67", "WHPREJUDICE67", "INCOME1F","UNION0" ,"CATHOLIC0",
"PROTESTANT0",   "EDYEARS6" , "FEMALE6","AGE6","LAGGOP","LAGNEITHER","WEAKDEM","INDDEM","INDEP","INDGOP"      
, "WEAKGOP","STRONGGOP" )

mf6$year <- "2012"
mf7$year <- "2012"
mf10$year <- "2016"
mf11$year <- "2016"
mf12$year <- "2016"

mf6$id <- rownames(mf6)
mf7$id <- rownames(mf7)
mf10$id <- rownames(mf10)
mf11$id <- rownames(mf11)
mf12$id <- rownames(mf12)

mfjoint <- rbind(mf6,mf7,mf10,mf11,mf12)

lout.all.wh.lm <- lmer(I(GOPVOTING==1) ~ WBPREJUDICE67*I(year=="2016")+WHPREJUDICE67*I(year=="2016")+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+LAGGOP+LAGNEITHER+WEAKDEM+INDDEM+INDEP+INDGOP+WEAKGOP+STRONGGOP+(1|id),data=mfjoint)
#### Appendix Table 13
texreg(lout.all.wh.lm,digits=3,stars=0.05)


###### Alternatives: Cruz, Rubio ----

res10rubio <- mnp(RUBIOVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=5000)
lout10rubio <- lm(RUBIOVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout10rubio)
colnames(mf)[1] <- "RUBIOVC10B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout10rubiob <- predict.mnp(res10rubio,newdata=mf1)
zero.10rubiob <- apply(apply(pout10rubiob$y,2,is.zero),2,mean)
quantile(zero.10rubiob,c(0.025,.5,.975))

one.10rubiob <- apply(apply(pout10rubiob$y,2,is.one),2,mean)
quantile(one.10rubiob,c(0.025,.5,.975))

two.10rubiob <- apply(apply(pout10rubiob$y,2,is.two),2,mean)
quantile(two.10rubiob,c(0.025,.5,.975))

pout10rubioc <- predict.mnp(res10rubio,newdata=mf2)
zero.10rubioc <- apply(apply(pout10rubioc$y,2,is.zero),2,mean)
quantile(zero.10rubioc,c(0.025,.5,.975))

one.10rubioc <- apply(apply(pout10rubioc$y,2,is.one),2,mean)
quantile(one.10rubioc,c(0.025,.5,.975))

two.10rubioc <- apply(apply(pout10rubioc$y,2,is.two),2,mean)
quantile(two.10rubioc,c(0.025,.5,.975))

gop.rubio.diff.10 <- two.10rubioc-two.10rubiob
quantile(gop.rubio.diff.10,c(0.025,.5,.975))

res10cruz <- mnp(CRUZVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw,n.draws=5000)
lout10cruz <- lm(CRUZVC10B ~ WBPREJUDICE67+WHPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+ROMNEYVO6A1+ROMNEYVO6A9+WEAKDEM6+INDDEM6+INDEP6+INDGOP6+WEAKGOP6+STRONGGOP6,data=dta.nhw)

mf <- model.frame(lout10cruz)
colnames(mf)[1] <- "CRUZVC10B"
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67,.8,na.rm=T)

pout10cruzb <- predict.mnp(res10cruz,newdata=mf1)
zero.10cruzb <- apply(apply(pout10cruzb$y,2,is.zero),2,mean)
quantile(zero.10cruzb,c(0.025,.5,.975))

one.10cruzb <- apply(apply(pout10cruzb$y,2,is.one),2,mean)
quantile(one.10cruzb,c(0.025,.5,.975))

two.10cruzb <- apply(apply(pout10cruzb$y,2,is.two),2,mean)
quantile(two.10cruzb,c(0.025,.5,.975))

pout10cruzc <- predict.mnp(res10cruz,newdata=mf2)
zero.10cruzc <- apply(apply(pout10cruzc$y,2,is.zero),2,mean)
quantile(zero.10cruzc,c(0.025,.5,.975))

one.10cruzc <- apply(apply(pout10cruzc$y,2,is.one),2,mean)
quantile(one.10cruzc,c(0.025,.5,.975))

two.10cruzc <- apply(apply(pout10cruzc$y,2,is.two),2,mean)
quantile(two.10cruzc,c(0.025,.5,.975))

gop.cruz.diff.10 <- two.10cruzc-two.10cruzb
quantile(gop.cruz.diff.10,c(0.025,.5,.975))

quantile(gop.diff.10.bp-gop.cruz.diff.10,c(0.025,.5,.975))

diff.cruz.trump.10 <- gop.diff.10-gop.cruz.diff.10
pval2 <- 1-sum(diff.cruz.trump.10 > 0)/length(diff.cruz.trump.10)

diff.rubio.trump.10 <- gop.diff.10-gop.rubio.diff.10
pval3 <- 1-sum(diff.rubio.trump.10 > 0)/length(diff.rubio.trump.10)

rmat10 <- matrix(NA,3,2)

# please note that 
rmat10[1,1] <- mean(gop.diff.10.bp)
rmat10[1,2] <- sd(gop.diff.10.bp)

rmat10[2,1] <- mean(gop.cruz.diff.10)
rmat10[2,2] <- sd(gop.cruz.diff.10)

rmat10[3,1] <- mean(gop.rubio.diff.10)
rmat10[3,2] <- sd(gop.rubio.diff.10)

pdf("trump-cruz-rubio-02132018-noarrows.pdf",width=8)

plot(x=1:3,y=rmat10[,1],pch=16,ylim=c(-.105,.125),xlim=c(0.5,3.5),cex.main=2,main="Change in GOP Support as \n Anti-Black Prejudice Rises",ylab="Shift in Prob. GOP Support",xlab="",cex.lab=1.3,xaxt="n")
for(i in 1:dim(rmat10)[1]){
	lines(x=c(i,i),y=c(rmat10[i,1]-rmat10[i,2]*1.96,rmat10[i,1]+rmat10[i,2]*1.96))	
}

text(x=1,y=.1175,paste("Trump \n vs. Clinton"))
brackets(x1=1,x2=2,y1=-0.05,y2=-0.05,h=-0.01)
text(x=1.5,y=-.065,paste("P=",round(pval2,digits=3),sep=""))

brackets(x1=1,x2=3,y1=-0.07,y2=-0.07,h=-0.01)
text(x=2,y=-.085,paste("P=",round(pval3,digits=3),sep=""))

text(x=2,y=.07,paste("Cruz \n vs. Clinton"))
text(x=3,y=.07,paste("Rubio \n vs. Clinton"))

abline(h=0)

dev.off()

#### primary 2016 ----
resp12 <- glm(TRUMP10P  ~ WBPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID6,data=dta.nhw,family=binomial(link=logit))

mf <- model.frame(resp12)
mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67[! dta.nhw$TRUMP10P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67[! dta.nhw$TRUMP10P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp12)$coef[,1],Sigma=vcov(resp12))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

trump.diff <- p2-p1
quantile(trump.diff,c(0.025,.5,.975))

resp12 <- glm(CLINTON10P  ~ WBPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID6,data=dta.nhw2,family=binomial(link=logit))

mf <- model.frame(resp12)
mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw2$WBPREJUDICE67[! dta.nhw2$SANDERS10P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw2$WBPREJUDICE67[! dta.nhw2$SANDERS10P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp12)$coef[,1],Sigma=vcov(resp12))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

clinton.diff <- p2-p1
quantile(clinton.diff,c(0.025,.5,.975))

resp12 <- glm(SANDERS10P  ~ WBPREJUDICE67+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID6,data=dta.nhw,family=binomial(link=logit))

mf <- model.frame(resp12)
mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67[! dta.nhw$SANDERS10P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE67 <- quantile(dta.nhw$WBPREJUDICE67[! dta.nhw$SANDERS10P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp12)$coef[,1],Sigma=vcov(resp12))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

sanders.diff <- p2-p1
quantile(sanders.diff,c(0.025,.5,.975))

#### primary 2008 ----

dta.nhw$GIULIANI1P <- 1*(dta.nhw$VP1A_REP_1=="Rudy Giuliani")
dta.nhw$THOMPSON1P <- 1*(dta.nhw$VP1A_REP_1=="Fred Thompson")
dta.nhw$ROMNEY1P <- 1*(dta.nhw$VP1A_REP_1=="Mitt Romney")
dta.nhw$MCCAIN1P <- 1*(dta.nhw$VP1A_REP_1=="John McCain")

dta.nhw2$OBAMA1P <- 1*(dta.nhw2$VP1A_DEM_1=="Barack Obama")

resp01 <- glm(ROMNEY1P  ~ WBPREJUDICE4+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID0,data=dta.nhw,family=binomial(link=logit))
mf <- model.frame(resp01)

mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$ROMNEY1P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$ROMNEY1P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp01)$coef[,1],Sigma=vcov(resp01))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

romney.diff <- p2-p1
quantile(romney.diff,c(0.025,.5,.975))

resp01 <- glm(MCCAIN1P  ~ WBPREJUDICE4+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID0,data=dta.nhw,family=binomial(link=logit))
mf <- model.frame(resp01)

mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$MCCAIN1P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$MCCAIN1P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp01)$coef[,1],Sigma=vcov(resp01))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

mccain.diff <- p2-p1
quantile(mccain.diff,c(0.025,.5,.975))

resp01 <- glm(THOMPSON1P  ~ WBPREJUDICE4+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID0,data=dta.nhw,family=binomial(link=logit))
mf <- model.frame(resp01)

mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$THOMPSON1P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$THOMPSON1P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp01)$coef[,1],Sigma=vcov(resp01))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

thompson.diff <- p2-p1
quantile(thompson.diff,c(0.025,.5,.975))

resp01 <- glm(GIULIANI1P  ~ WBPREJUDICE4+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID0,data=dta.nhw,family=binomial(link=logit))
mf <- model.frame(resp01)

mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$GIULIANI1P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE4 <- quantile(dta.nhw$WBPREJUDICE4[! dta.nhw$GIULIANI1P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp01)$coef[,1],Sigma=vcov(resp01))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

giuliani.diff <- p2-p1
quantile(giuliani.diff,c(0.025,.5,.975))

### OBAMA

resp01 <- glm(OBAMA1P  ~ WBPREJUDICE4+INCOME1F+UNION0+CATHOLIC0+PROTESTANT0+EDYEARS6+FEMALE6+AGE6+PID0,data=dta.nhw2,family=binomial(link=logit))
mf <- model.frame(resp01)

mf[,1] <- rep(1,dim(mf)[1])
mf1 <- mf2 <- mf
mf1$WBPREJUDICE4 <- quantile(dta.nhw2$WBPREJUDICE4[! dta.nhw2$OBAMA1P %in% c(NA)],.2,na.rm=T)
mf2$WBPREJUDICE4 <- quantile(dta.nhw2$WBPREJUDICE4[! dta.nhw2$OBAMA1P %in% c(NA)],.8,na.rm=T)

mf1.vec <- apply(mf1,2,FUN=median)
mf2.vec <- apply(mf2,2,FUN=median)

M <- 10000
coefs <- mvrnorm(M,mu=summary(resp01)$coef[,1],Sigma=vcov(resp01))

p1 <- 1/(1+exp(-coefs%*%mf1.vec))
p2 <- 1/(1+exp(-coefs%*%mf2.vec))

obama.diff <- p2-p1
quantile(obama.diff,c(0.025,.5,.975))

rmat20 <- matrix(NA,8,2)

rmat20[1,1] <- mean(trump.diff)
rmat20[1,2] <- sd(trump.diff)

rmat20[2,1] <- mean(thompson.diff)
rmat20[2,2] <- sd(thompson.diff)

rmat20[3,1] <- mean(giuliani.diff)
rmat20[3,2] <- sd(giuliani.diff)

rmat20[4,1] <- mean(mccain.diff)
rmat20[4,2] <- sd(mccain.diff)

rmat20[5,1] <- mean(romney.diff)
rmat20[5,2] <- sd(romney.diff)

rmat20[6,1] <- mean(clinton.diff)
rmat20[6,2] <- sd(clinton.diff)

rmat20[7,1] <- mean(sanders.diff)
rmat20[7,2] <- sd(sanders.diff)

rmat20[8,1] <- mean(obama.diff)
rmat20[8,2] <- sd(obama.diff)

#pdf("primaries-07312019-noarrows.pdf",width=12)
#jpeg("primaries-07312019-noarrows.jpg",width=720)

plot(x=1:8,y=rmat20[,1],pch=16,ylim=c(-.125,.205),xlim=c(0.5,8.5),cex.main=2,main="Change in Primary Support as \n Anti-Black Prejudice Rises",ylab="Shift in Prob. Support",xlab="",cex.lab=1.3,xaxt="n")
for(i in 1:dim(rmat20)[1]){
	lines(x=c(i,i),y=c(rmat20[i,1]-rmat20[i,2]*1.96,rmat20[i,1]+rmat20[i,2]*1.96))	
}

text(x=1,y=.18,paste("Trump \n 2016 Primary"))
text(x=2,y=.18,paste("Thompson \n 2008 Primary"))
text(x=3,y=.18,paste("Giuliani \n 2008 Primary"))
text(x=4,y=.18,paste("McCain \n 2008 Primary"))
text(x=5,y=.18,paste("Romney \n 2008 Primary"))
text(x=6,y=.18,paste("Clinton \n 2008 Primary"))
text(x=7,y=.18,paste("Sanders \n 2016 Primary"))
text(x=8,y=.18,paste("Obama \n 2008 Primary"))

abline(h=0)

dev.off()

